home *** CD-ROM | disk | FTP | other *** search
/ OpenGL Superbible (2nd Edition) / OpenGL SuperBible e2.iso / tools / GLUT-3.7 / PROGS / ADVANCED / IZOOM.C < prev    next >
Encoding:
C/C++ Source or Header  |  1998-08-12  |  15.0 KB  |  672 lines
  1. /**
  2.  **    izoom- 
  3.  **        Magnify or minify a picture with or without filtering.  The 
  4.  **    filtered method is one pass, uses 2-d convolution, and is optimized 
  5.  **    by integer arithmetic and precomputation of filter coeffs.
  6.  **
  7.  **                    Paul Haeberli - 1988
  8.  **/
  9. #include <stdlib.h>
  10. #include <stdio.h>
  11. #include <string.h>
  12. #include <math.h>
  13. #include "izoom.h"
  14.  
  15. #ifdef _WIN32
  16. #pragma warning (disable:4244)          /* disable bogus conversion warnings */
  17. #endif
  18.  
  19. typedef struct filtinteg {
  20.   float rad, min, max;
  21.   float *tab;
  22. } filtinteg;
  23.  
  24. float flerp(float f0, float f1, float p);
  25.  
  26. #define GRIDTOFLOAT(pos,n)    (((pos)+0.5)/(n))
  27. #define FLOATTOGRID(pos,n)    ((pos)*(n))
  28. #define SHIFT             12
  29. #define ONE             (1<<SHIFT)
  30. #define EPSILON            0.0001
  31. #define FILTERRAD        (blurfactor*shape->rad)
  32. #define FILTTABSIZE        250
  33.  
  34. static void makexmap(short *abuf, short *xmap[], int anx, int bnx);
  35. static void setintrow(int *buf, int val, int n);
  36. static void xscalebuf(short *xmap[], short *bbuf, int bnx);
  37. static void addrow(int *iptr, short *sptr, int w, int n);
  38. static void divrow(int *iptr, short *sptr, int tot, int n);
  39. static FILTER *makefilt(short *abuf, int anx, int bnx, int *maxn);
  40. static void freefilt(FILTER * filt, int n);
  41. static void applyxfilt(short *bbuf, FILTER * xfilt, int bnx);
  42. float filterinteg(float bmin, float bmax, float blurf);
  43. static void mitchellinit(float b, float c);
  44. static void clamprow(short *iptr, short *optr, int n);
  45.  
  46. float filt_box(float x);
  47. float filt_triangle(float x);
  48. float filt_quadratic(float x);
  49. float filt_mitchell(float x);
  50. float filt_gaussian(float x);
  51.  
  52. static int (*xfiltfunc) (short *, int);
  53. static float blurfactor;
  54. int izoomdebug;
  55.  
  56. static filtinteg *shapeBOX;
  57. static filtinteg *shapeTRIANGLE;
  58. static filtinteg *shapeQUADRATIC;
  59. static filtinteg *shapeMITCHELL;
  60. static filtinteg *shapeGAUSSIAN;
  61. static filtinteg *shape;
  62.  
  63. static filtinteg *
  64. integrate(float (*filtfunc) (float), float rad)
  65. {
  66.   int i;
  67.   float del, x, min, max;
  68.   double tot;
  69.   filtinteg *filt;
  70.  
  71.   min = -rad;
  72.   max = rad;
  73.   del = 2 * rad;
  74.   tot = 0.0;
  75.   filt = (filtinteg *) malloc(sizeof(filtinteg));
  76.   filt->rad = rad;
  77.   filt->min = min;
  78.   filt->max = max;
  79.   filt->tab = (float *) malloc(FILTTABSIZE * sizeof(float));
  80.   for (i = 0; i < FILTTABSIZE; i++) {
  81.     x = min + (del * i / (FILTTABSIZE - 1.0));
  82.     tot = tot + filtfunc(x);
  83.     filt->tab[i] = tot;
  84.   }
  85.   for (i = 0; i < FILTTABSIZE; i++)
  86.     filt->tab[i] /= tot;
  87.   return filt;
  88. }
  89.  
  90. float 
  91. filterinteg(float bmin, float bmax, float blurf)
  92. {
  93.   int i1, i2;
  94.   float f1, f2;
  95.   float *tab;
  96.   float mult;
  97.  
  98.   bmin /= blurf;
  99.   bmax /= blurf;
  100.   tab = shape->tab;
  101.   mult = (FILTTABSIZE - 1.0) / (2.0 * shape->rad);
  102.  
  103.   f1 = ((bmin - shape->min) * mult);
  104.   i1 = floor(f1);
  105.   f1 = f1 - i1;
  106.   if (i1 < 0)
  107.     f1 = 0.0;
  108.   else if (i1 >= (FILTTABSIZE - 1))
  109.     f1 = 1.0;
  110.   else
  111.     f1 = flerp(tab[i1], tab[i1 + 1], f1);
  112.  
  113.   f2 = ((bmax - shape->min) * mult);
  114.   i2 = floor(f2);
  115.   f2 = f2 - i2;
  116.   if (i2 < 0)
  117.     f2 = 0.0;
  118.   else if (i2 >= (FILTTABSIZE - 1))
  119.     f2 = 1.0;
  120.   else
  121.     f2 = flerp(tab[i2], tab[i2 + 1], f2);
  122.   return f2 - f1;
  123. }
  124.  
  125. void
  126. setfiltertype(int filttype)
  127. {
  128.   switch (filttype) {
  129.   case IMPULSE:
  130.     shape = 0;
  131.     break;
  132.   case BOX:
  133.     if (!shapeBOX)
  134.       shapeBOX = integrate(filt_box, 0.5);
  135.     shape = shapeBOX;
  136.     break;
  137.   case TRIANGLE:
  138.     if (!shapeTRIANGLE)
  139.       shapeTRIANGLE = integrate(filt_triangle, 1.0);
  140.     shape = shapeTRIANGLE;
  141.     break;
  142.   case QUADRATIC:
  143.     if (!shapeQUADRATIC)
  144.       shapeQUADRATIC = integrate(filt_quadratic, 1.5);
  145.     shape = shapeQUADRATIC;
  146.     break;
  147.   case MITCHELL:
  148.     if (!shapeMITCHELL)
  149.       shapeMITCHELL = integrate(filt_mitchell, 2.0);
  150.     shape = shapeMITCHELL;
  151.     break;
  152.   case GAUSSIAN:
  153.     if (!shapeGAUSSIAN)
  154.       shapeGAUSSIAN = integrate(filt_gaussian, 1.5);
  155.     shape = shapeGAUSSIAN;
  156.     break;
  157.   }
  158. }
  159.  
  160. void
  161. copyimage(getfunc_t getfunc, getfunc_t putfunc, int nx, int ny)
  162. {
  163.   int y;
  164.   short *abuf;
  165.  
  166.   abuf = (short *) malloc(nx * sizeof(short));
  167.   for (y = 0; y < ny; y++) {
  168.     getfunc(abuf, y);
  169.     putfunc(abuf, y);
  170.   }
  171.   free(abuf);
  172. }
  173.  
  174. /* general zoom follows */
  175. zoom *
  176. newzoom(getfunc_t getfunc, int anx, int any, int bnx, int bny, int filttype, float blur)
  177. {
  178.   zoom *z;
  179.   int i;
  180.  
  181.   setfiltertype(filttype);
  182.   z = (zoom *) malloc(sizeof(zoom));
  183.   z->getfunc = getfunc;
  184.   z->abuf = (short *) malloc(anx * sizeof(short));
  185.   z->bbuf = (short *) malloc(bnx * sizeof(short));
  186.   z->anx = anx;
  187.   z->any = any;
  188.   z->bnx = bnx;
  189.   z->bny = bny;
  190.   z->curay = -1;
  191.   z->y = 0;
  192.   z->type = filttype;
  193.   if (filttype == IMPULSE) {
  194.     if (z->anx != z->bnx) {
  195.       z->xmap = (short **) malloc(z->bnx * sizeof(short *));
  196.       makexmap(z->abuf, z->xmap, z->anx, z->bnx);
  197.     }
  198.   } else {
  199.     blurfactor = blur;
  200.     if (filttype == MITCHELL)
  201.       z->clamp = 1;
  202.     else
  203.       z->clamp = 0;
  204.     z->tbuf = (short *) malloc(bnx * sizeof(short));
  205.     z->xfilt = makefilt(z->abuf, anx, bnx, &z->nrows);
  206.     z->yfilt = makefilt(0, any, bny, &z->nrows);
  207.     z->filtrows = (short **) malloc(z->nrows * sizeof(short *));
  208.     for (i = 0; i < z->nrows; i++)
  209.       z->filtrows[i] = (short *) malloc(z->bnx * sizeof(short));
  210.     z->accrow = (int *) malloc(z->bnx * sizeof(int));
  211.     z->ay = 0;
  212.   }
  213.   return z;
  214. }
  215.  
  216. void
  217. getzoomrow(zoom * z, short *buf, int y)
  218. {
  219.   float fy;
  220.   int ay;
  221.   FILTER *f;
  222.   int i, max;
  223.   short *row;
  224.  
  225.   if (y == 0) {
  226.     z->curay = -1;
  227.     z->y = 0;
  228.     z->ay = 0;
  229.   }
  230.   if (z->type == IMPULSE) {
  231.     fy = GRIDTOFLOAT(z->y, z->bny);
  232.     ay = FLOATTOGRID(fy, z->any);
  233.     if (z->anx == z->bnx) {
  234.       if (z->curay != ay) {
  235.         z->getfunc(z->abuf, ay);
  236.         z->curay = ay;
  237.         if (xfiltfunc)
  238.           xfiltfunc(z->abuf, z->bnx);
  239.       }
  240.       memcpy(buf, z->abuf, z->bnx * sizeof(short));
  241.     } else {
  242.       if (z->curay != ay) {
  243.         z->getfunc(z->abuf, ay);
  244.         xscalebuf(z->xmap, z->bbuf, z->bnx);
  245.         z->curay = ay;
  246.         if (xfiltfunc)
  247.           xfiltfunc(z->bbuf, z->bnx);
  248.       }
  249.       memcpy(buf, z->bbuf, z->bnx * sizeof(short));
  250.     }
  251.   } else if (z->any == 1 && z->bny == 1) {
  252.     z->getfunc(z->abuf, z->ay++);
  253.     applyxfilt(z->filtrows[0], z->xfilt, z->bnx);
  254.     if (xfiltfunc)
  255.       xfiltfunc(z->filtrows[0], z->bnx);
  256.     if (z->clamp) {
  257.       clamprow(z->filtrows[0], z->tbuf, z->bnx);
  258.       memcpy(buf, z->tbuf, z->bnx * sizeof(short));
  259.     } else {
  260.       memcpy(buf, z->filtrows[0], z->bnx * sizeof(short));
  261.     }
  262.   } else {
  263.     f = z->yfilt + z->y;
  264.     max = (int) (sizeof(f->dat) / sizeof(short) + (f->n - 1));
  265.     while (z->ay <= max) {
  266.       z->getfunc(z->abuf, z->ay++);
  267.       row = z->filtrows[0];
  268.       for (i = 0; i < (z->nrows - 1); i++)
  269.         z->filtrows[i] = z->filtrows[i + 1];
  270.       z->filtrows[z->nrows - 1] = row;
  271.       applyxfilt(z->filtrows[z->nrows - 1], z->xfilt, z->bnx);
  272.       if (xfiltfunc)
  273.         xfiltfunc(z->filtrows[z->nrows - 1], z->bnx);
  274.     }
  275.     if (f->n == 1) {
  276.       if (z->clamp) {
  277.         clamprow(z->filtrows[z->nrows - 1], z->tbuf, z->bnx);
  278.         memcpy(buf, z->tbuf, z->bnx * sizeof(short));
  279.       } else {
  280.         memcpy(buf, z->filtrows[z->nrows - 1], z->bnx * sizeof(short));
  281.       }
  282.     } else {
  283.       setintrow(z->accrow, f->halftotw, z->bnx);
  284.       for (i = 0; i < f->n; i++)
  285.         addrow(z->accrow, z->filtrows[i + (z->nrows - 1) - (f->n - 1)],
  286.           f->w[i], z->bnx);
  287.       divrow(z->accrow, z->bbuf, f->totw, z->bnx);
  288.       if (z->clamp) {
  289.         clamprow(z->bbuf, z->tbuf, z->bnx);
  290.         memcpy(buf, z->tbuf, z->bnx * sizeof(short));
  291.       } else {
  292.         memcpy(buf, z->bbuf, z->bnx * sizeof(short));
  293.       }
  294.     }
  295.   }
  296.   z->y++;
  297. }
  298.  
  299. static void
  300. setintrow(int *buf, int val, int n)
  301. {
  302.   while (n >= 8) {
  303.     buf[0] = val;
  304.     buf[1] = val;
  305.     buf[2] = val;
  306.     buf[3] = val;
  307.     buf[4] = val;
  308.     buf[5] = val;
  309.     buf[6] = val;
  310.     buf[7] = val;
  311.     buf += 8;
  312.     n -= 8;
  313.   }
  314.   while (n--)
  315.     *buf++ = val;
  316. }
  317.  
  318. void
  319. freezoom(zoom * z)
  320. {
  321.   int i;
  322.  
  323.   if (z->type == IMPULSE) {
  324.     if (z->anx != z->bnx)
  325.       free(z->xmap);
  326.   } else {
  327.     freefilt(z->xfilt, z->bnx);
  328.     freefilt(z->yfilt, z->bny);
  329.     free(z->tbuf);
  330.     for (i = 0; i < z->nrows; i++)
  331.       free(z->filtrows[i]);
  332.     free(z->filtrows);
  333.     free(z->accrow);
  334.   }
  335.   free(z->abuf);
  336.   free(z->bbuf);
  337.   free(z);
  338.  
  339. }
  340.  
  341. void
  342. filterzoom(getfunc_t getfunc, getfunc_t putfunc, int anx, int any, int bnx, int bny, int filttype, float blur)
  343. {
  344.   zoom *z;
  345.   int y;
  346.   short *buf;
  347.  
  348.   buf = (short *) malloc(bnx * sizeof(short));
  349.   z = newzoom(getfunc, anx, any, bnx, bny, filttype, blur);
  350.   for (y = 0; y < bny; y++) {
  351.     getzoomrow(z, buf, y);
  352.     putfunc(buf, y);
  353.   }
  354.   freezoom(z);
  355.   free(buf);
  356. }
  357.  
  358. /* impulse zoom utilities */
  359. static void
  360. makexmap(short *abuf, short *xmap[], int anx, int bnx)
  361. {
  362.   int x, ax;
  363.   float fx;
  364.  
  365.   for (x = 0; x < bnx; x++) {
  366.     fx = GRIDTOFLOAT(x, bnx);
  367.     ax = FLOATTOGRID(fx, anx);
  368.     xmap[x] = abuf + ax;
  369.   }
  370. }
  371.  
  372. static void
  373. xscalebuf(short *xmap[], short *bbuf, int bnx)
  374. {
  375.   while (bnx >= 8) {
  376.     bbuf[0] = *(xmap[0]);
  377.     bbuf[1] = *(xmap[1]);
  378.     bbuf[2] = *(xmap[2]);
  379.     bbuf[3] = *(xmap[3]);
  380.     bbuf[4] = *(xmap[4]);
  381.     bbuf[5] = *(xmap[5]);
  382.     bbuf[6] = *(xmap[6]);
  383.     bbuf[7] = *(xmap[7]);
  384.     bbuf += 8;
  385.     xmap += 8;
  386.     bnx -= 8;
  387.   }
  388.   while (bnx--)
  389.     *bbuf++ = *(*xmap++);
  390. }
  391.  
  392. void
  393. zoomxfilt(int (*filtfunc) (short *, int))
  394. {
  395.   xfiltfunc = filtfunc;
  396. }
  397.  
  398. /* filter zoom utilities */
  399. static void
  400. addrow(int *iptr, short *sptr, int w, int n)
  401. {
  402.   while (n >= 8) {
  403.     iptr[0] += (w * sptr[0]);
  404.     iptr[1] += (w * sptr[1]);
  405.     iptr[2] += (w * sptr[2]);
  406.     iptr[3] += (w * sptr[3]);
  407.     iptr[4] += (w * sptr[4]);
  408.     iptr[5] += (w * sptr[5]);
  409.     iptr[6] += (w * sptr[6]);
  410.     iptr[7] += (w * sptr[7]);
  411.     iptr += 8;
  412.     sptr += 8;
  413.     n -= 8;
  414.   }
  415.   while (n--)
  416.     *iptr++ += (w * *sptr++);
  417. }
  418.  
  419. static void
  420. divrow(int *iptr, short *sptr, int tot, int n)
  421. {
  422.   while (n >= 8) {
  423.     sptr[0] = iptr[0] / tot;
  424.     sptr[1] = iptr[1] / tot;
  425.     sptr[2] = iptr[2] / tot;
  426.     sptr[3] = iptr[3] / tot;
  427.     sptr[4] = iptr[4] / tot;
  428.     sptr[5] = iptr[5] / tot;
  429.     sptr[6] = iptr[6] / tot;
  430.     sptr[7] = iptr[7] / tot;
  431.     sptr += 8;
  432.     iptr += 8;
  433.     n -= 8;
  434.   }
  435.   while (n--)
  436.     *sptr++ = (*iptr++) / tot;
  437. }
  438.  
  439. static FILTER *
  440. makefilt(short *abuf, int anx, int bnx, int *maxn)
  441. {
  442.   FILTER *f, *filter;
  443.   int x, n;
  444.   float bmin, bmax, bcent, brad;
  445.   float fmin, fmax, acent, arad;
  446.   int amin, amax;
  447.   float coverscale;
  448.  
  449.   if (izoomdebug)
  450.     fprintf(stderr, "makefilt\n");
  451.   f = filter = (FILTER *) malloc(bnx * sizeof(FILTER));
  452.   *maxn = 0;
  453.   if (bnx < anx) {
  454.     coverscale = ((float) anx / bnx * ONE) / 2.0;
  455.     brad = FILTERRAD / bnx;
  456.     for (x = 0; x < bnx; x++) {
  457.       bcent = ((float) x + 0.5) / bnx;
  458.       amin = floor((bcent - brad) * anx + EPSILON);
  459.       amax = floor((bcent + brad) * anx - EPSILON);
  460.       if (amin < 0)
  461.         amin = 0;
  462.       if (amax >= anx)
  463.         amax = anx - 1;
  464.       f->n = 1 + amax - amin;
  465.       f->dat = abuf + amin;
  466.       f->w = (short *) malloc(f->n * sizeof(short));
  467.       f->totw = 0;
  468.       if (izoomdebug)
  469.         fprintf(stderr, "| ");
  470.       for (n = 0; n < f->n; n++) {
  471.         bmin = bnx * ((((float) amin + n) / anx) - bcent);
  472.         bmax = bnx * ((((float) amin + n + 1) / anx) - bcent);
  473.         f->w[n] = floor((coverscale * filterinteg(bmin, bmax, blurfactor)) + 0.5);
  474.         if (izoomdebug)
  475.           fprintf(stderr, "%d ", f->w[n]);
  476.         f->totw += f->w[n];
  477.       }
  478.       f->halftotw = f->totw / 2;
  479.       if (f->n > *maxn)
  480.         *maxn = f->n;
  481.       f++;
  482.     }
  483.   } else {
  484.     coverscale = ((float) bnx / anx * ONE) / 2.0;
  485.     arad = FILTERRAD / anx;
  486.     for (x = 0; x < bnx; x++) {
  487.       bmin = ((float) x) / bnx;
  488.       bmax = ((float) x + 1.0) / bnx;
  489.       amin = floor((bmin - arad) * anx + (0.5 + EPSILON));
  490.       amax = floor((bmax + arad) * anx - (0.5 + EPSILON));
  491.       if (amin < 0)
  492.         amin = 0;
  493.       if (amax >= anx)
  494.         amax = anx - 1;
  495.       f->n = 1 + amax - amin;
  496.       f->dat = abuf + amin;
  497.       f->w = (short *) malloc(f->n * sizeof(short));
  498.       f->totw = 0;
  499.       if (izoomdebug)
  500.         fprintf(stderr, "| ");
  501.       for (n = 0; n < f->n; n++) {
  502.         acent = (amin + n + 0.5) / anx;
  503.         fmin = anx * (bmin - acent);
  504.         fmax = anx * (bmax - acent);
  505.         f->w[n] = floor((coverscale * filterinteg(fmin, fmax, blurfactor)) + 0.5);
  506.         if (izoomdebug)
  507.           fprintf(stderr, "%d ", f->w[n]);
  508.         f->totw += f->w[n];
  509.       }
  510.       f->halftotw = f->totw / 2;
  511.       if (f->n > *maxn)
  512.         *maxn = f->n;
  513.       f++;
  514.     }
  515.   }
  516.   if (izoomdebug)
  517.     fprintf(stderr, "|\n");
  518.   return filter;
  519. }
  520.  
  521. static void
  522. freefilt(FILTER * filt, int n)
  523. {
  524.   FILTER *f;
  525.  
  526.   f = filt;
  527.   while (n--) {
  528.     free(f->w);
  529.     f++;
  530.   }
  531.   free(filt);
  532. }
  533.  
  534. static void
  535. applyxfilt(short *bbuf, FILTER * xfilt, int bnx)
  536. {
  537.   short *w;
  538.   short *dptr;
  539.   int n, val;
  540.  
  541.   while (bnx--) {
  542.     if ((n = xfilt->n) == 1) {
  543.       *bbuf++ = *xfilt->dat;
  544.     } else {
  545.       w = xfilt->w;
  546.       dptr = xfilt->dat;
  547.       val = xfilt->halftotw;
  548.       n = xfilt->n;
  549.       while (n--)
  550.         val += *w++ * *dptr++;
  551.       *bbuf++ = val / xfilt->totw;
  552.     }
  553.     xfilt++;
  554.   }
  555. }
  556.  
  557. /* filter shape functions follow */
  558. float 
  559. filt_box(float x)
  560. {
  561.   if (x < -0.5)
  562.     return 0.0;
  563.   if (x < 0.5)
  564.     return 1.0;
  565.   return 0.0;
  566. }
  567.  
  568. float 
  569. filt_triangle(float x)
  570. {
  571.   if (x < -1.0)
  572.     return 0.0;
  573.   if (x < 0.0)
  574.     return 1.0 + x;
  575.   if (x < 1.0)
  576.     return 1.0 - x;
  577.   return 0.0;
  578. }
  579.  
  580. float 
  581. filt_quadratic(float x)
  582. {
  583.   if (x < -1.5)
  584.     return 0.0;
  585.   if (x < -0.5)
  586.     return 0.5 * (x + 1.5) * (x + 1.5);
  587.   if (x < 0.5)
  588.     return 0.75 - (x * x);
  589.   if (x < 1.5)
  590.     return 0.5 * (x - 1.5) * (x - 1.5);
  591.   return 0.0;
  592. }
  593.  
  594. static float p0, p2, p3, q0, q1, q2, q3;
  595.  
  596. /* see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  597.    SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.  */
  598.  
  599. float 
  600. filt_mitchell(float x)
  601. {                       /* Mitchell & Netravali's two-param cubic */
  602.   static int mitfirsted;
  603.  
  604.   if (!mitfirsted) {
  605.     mitchellinit(1.0f / 3.0f, 1.0f / 3.0f);
  606.     mitfirsted = 1;
  607.   }
  608.   if (x < -2.0)
  609.     return 0.0;
  610.   if (x < -1.0)
  611.     return (q0 - x * (q1 - x * (q2 - x * q3)));
  612.   if (x < 0.0)
  613.     return (p0 + x * x * (p2 - x * p3));
  614.   if (x < 1.0)
  615.     return (p0 + x * x * (p2 + x * p3));
  616.   if (x < 2.0)
  617.     return (q0 + x * (q1 + x * (q2 + x * q3)));
  618.   return 0.0;
  619. }
  620.  
  621. static void
  622. mitchellinit(float b, float c)
  623. {
  624.   p0 = (6.0 - 2.0 * b) / 6.0;
  625.   p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0;
  626.   p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0;
  627.   q0 = (8.0 * b + 24.0 * c) / 6.0;
  628.   q1 = (-12.0 * b - 48.0 * c) / 6.0;
  629.   q2 = (6.0 * b + 30.0 * c) / 6.0;
  630.   q3 = (-b - 6.0 * c) / 6.0;
  631. }
  632.  
  633. #define NARROWNESS    1.5
  634.  
  635. float 
  636. filt_gaussian(float x)
  637. {
  638.   x = x * NARROWNESS;
  639.   return (1.0 / exp(x * x) - 1.0 / exp(1.5 * NARROWNESS * 1.5 * NARROWNESS));
  640. }
  641.  
  642. float 
  643. flerp(float f0, float f1, float p)
  644. {
  645.   return ((f0 * (1.0 - p)) + (f1 * p));
  646. }
  647.  
  648. #define DOCLAMP(iptr,optr)    *(optr) = ((*(iptr)<0) ? 0 : (*(iptr)>255) ? 255 : *(iptr))
  649.  
  650. static void
  651. clamprow(short *iptr, short *optr, int n)
  652. {
  653.   while (n >= 8) {
  654.     DOCLAMP(iptr + 0, optr + 0);
  655.     DOCLAMP(iptr + 1, optr + 1);
  656.     DOCLAMP(iptr + 2, optr + 2);
  657.     DOCLAMP(iptr + 3, optr + 3);
  658.     DOCLAMP(iptr + 4, optr + 4);
  659.     DOCLAMP(iptr + 5, optr + 5);
  660.     DOCLAMP(iptr + 6, optr + 6);
  661.     DOCLAMP(iptr + 7, optr + 7);
  662.     iptr += 8;
  663.     optr += 8;
  664.     n -= 8;
  665.   }
  666.   while (n--) {
  667.     DOCLAMP(iptr, optr);
  668.     iptr++;
  669.     optr++;
  670.   }
  671. }
  672.